#ifndef VOLUME_DATA_H
#define VOLUME_DATA_H

#include <G3D/G3DAll.h>
#include <GLG3D/GLG3D.h>
#include <fstream>

#include "MarchingCubes.h"

template<typename T>
struct VolumeData
{
public:
	T minval, maxval;
	Vector3int32 dim;
	Matrix4 v2w;
	Matrix4 w2v;
	Array<T> data;
	Array<Vector3> vertexArray;
	Array<Vector3> normalArray;
	Array<Color4uint8> colorArray;
	VertexBufferRef vbo;
	VertexRange vb;
	VertexRange nb;
	VertexRange cb;
	TextureRef vtex;

public:
	void load(const std::string &file, int width, int height, int depth)
	{
		std::string txt = file + ".txt";
		{
			std::ifstream fin(txt.c_str(), std::ios::in | std::ios::binary);
			if(!fin.is_open())
				return;
			std::string buf;
			fin >> buf;
			Vector4 row;
			fin >> row.x >> row.y >> row.z >> row.w;
			v2w.setRow(0, row);
			fin >> row.x >> row.y >> row.z >> row.w;
			v2w.setRow(1, row);
			fin >> row.x >> row.y >> row.z >> row.w;
			v2w.setRow(2, row);
			fin >> row.x >> row.y >> row.z >> row.w;
			v2w.setRow(3, row);
			w2v = v2w.inverse();
		}

		std::string raw = file + ".raw";
		{
			std::ifstream fin(raw.c_str(), std::ios::in | std::ios::binary);
			if(!fin.is_open())
				return;
			int size = width * height * depth;
			dim.x = width;
			dim.y = height;
			dim.z = depth;
			data.resize(size);
			fin.read((char *)&data[0], size * sizeof(T));
			fin.close();
		}

		findMinMax();
	}

	void save(const std::string &file)
	{
		std::ofstream fout(file.c_str(), std::ios::out | std::ios::binary);
		if(!fout.is_open())
			return;
		fout.write((char *)&data[0], data.size() * sizeof(T));
		fout.close();
	}

	void findMinMax()
	{
		if(data.size() == 0)
		{
			minval = T();
			maxval = T();
			return;
		}

		minval = data[0];
		maxval = data[0];
		for(int i = 1; i < data.size(); ++i)
		{
			minval = min(minval, data[i]);
			maxval = max(maxval, data[i]);
		}
	}

	void smoothScalar(float mu)
	{
		Array<T> temp(data);
		int nx = dim.x, nxy = dim.x * dim.y;
		for(int it = 0; it < 16; ++it)
		{
			for(int k = 1; k < dim.z - 1; ++k) for(int j = 1; j < dim.y - 1; ++j) for(int i = 1; i < dim.x - 1; ++i)
			{
				int idx = i + j * nx + k * nxy;
				data[idx] = T((temp[idx] + (T(0) + data[idx - 1] + data[idx + 1] + data[idx - nx] + data[idx + nx] + data[idx - nxy] + data[idx + nxy]) * mu) / (1.0f + 6.0f * mu));
			}
		}
	}

	void smoothVector(float mu)
	{
		Array<T> temp(data);
		int nx = dim.x, nxy = dim.x * dim.y;
		for(int it = 0; it < 16; ++it)
		{
			for(int k = 1; k < dim.z - 1; ++k) for(int j = 1; j < dim.y - 1; ++j) for(int i = 1; i < dim.x - 1; ++i)
			{
				int idx = i + j * nx + k * nxy;
				data[idx] = T((temp[idx] + (T() + data[idx - 1] + data[idx + 1] + data[idx - nx] + data[idx + nx] + data[idx - nxy] + data[idx + nxy]) * mu) / (1.0f + 6.0f * mu));
			}
		}
	}

	void createIsoSurface(T isoValue)
	{
		int offset = vertexArray.size();
		MarchingCubes<T>::generateMesh(dim.x, dim.y, dim.z, data, isoValue, vertexArray, normalArray);
		for(int i = offset; i < vertexArray.size(); ++i)
		{
			vertexArray[i] = v2w.homoMul(vertexArray[i], 1.0f);
			normalArray[i] = (v2w.upper3x3() * normalArray[i]).directionOrZero();
		}

	}

	void updateVertexBuffer()
	{
		vbo = VertexBuffer::create((vertexArray.size() + normalArray.size() + colorArray.size()) * sizeof(Vector3) + 32, VertexBuffer::WRITE_ONCE, VertexBuffer::DATA);
		if(vertexArray.size() > 0)
			vb = VertexRange(vertexArray, vbo);
		if(normalArray.size() > 0)
			nb = VertexRange(normalArray, vbo);
		if(colorArray.size() > 0)
			cb = VertexRange(colorArray, vbo);
	}

	void updateVolumeTexture();

	void renderIsoSurface(RenderDevice *rd)
	{
		rd->pushState();
		//rd->setRenderMode(RenderDevice::RENDER_WIREFRAME);
		//rd->setBlendFunc(RenderDevice::BLEND_SRC_ALPHA, RenderDevice::BLEND_ONE_MINUS_SRC_ALPHA);
		rd->beginIndexedPrimitives();
		rd->setVertexArray(vb);
		if(normalArray.size() > 0)
			rd->setNormalArray(nb);
		if(colorArray.size() > 0)
			rd->setColorArray(cb);
		rd->sendSequentialIndices(RenderDevice::Primitive::TRIANGLES, vertexArray.size());
		rd->endIndexedPrimitives();
		rd->popState();
	}

	void renderVectorField(RenderDevice *rd)
	{
		rd->pushState();
		glDisable(GL_LINE_SMOOTH);
		//rd->setRenderMode(RenderDevice::RENDER_WIREFRAME);
		//rd->setBlendFunc(RenderDevice::BLEND_SRC_ALPHA, RenderDevice::BLEND_ONE_MINUS_SRC_ALPHA);
		//rd->setDepthTest(RenderDevice::DEPTH_ALWAYS_PASS);
		//rd->setDepthWrite(false);
		rd->setLineWidth(1.0f);
		rd->setShadeMode(RenderDevice::SHADE_SMOOTH);
		rd->beginIndexedPrimitives();
		rd->setVertexArray(vb);
		if(normalArray.size() > 0)
			rd->setNormalArray(nb);
		if(colorArray.size() > 0)
			rd->setColorArray(cb);
		rd->sendSequentialIndices(RenderDevice::Primitive::LINES, vertexArray.size());
		rd->endIndexedPrimitives();
		glEnable(GL_LINE_SMOOTH);
		rd->popState();
	}

	inline T sample(float x, float y, float z)
	{
		x = clamp(x, 0, dim.x - 1);
		y = clamp(y, 0, dim.y - 1);
		z = clamp(z, 0, dim.z - 1);
		int xi = floor(x);
		int yi = floor(y);
		int zi = floor(z);
		float xf = x - xi;
		float yf = y - yi;
		float zf = z - zi;
		if(xi >= dim.x - 1) { --xi; xf += 1.0f; }
		if(yi >= dim.y - 1) { --yi; yf += 1.0f; }
		if(zi >= dim.z - 1) { --zi; zf += 1.0f; }
		int i0 = xi + yi * dim.x + zi * dim.x * dim.y;
		int i1 = i0 + 1;
		int i2 = i0 + dim.x;
		int i3 = i2 + 1;
		int i4 = i0 + dim.x * dim.y;
		int i5 = i4 + 1;
		int i6 = i4 + dim.x;
		int i7 = i6 + 1;
		T v01 = data[i0] * (1.0f - xf) + data[i1] * xf;
		T v23 = data[i2] * (1.0f - xf) + data[i3] * xf;
		T v45 = data[i4] * (1.0f - xf) + data[i5] * xf;
		T v67 = data[i6] * (1.0f - xf) + data[i7] * xf;
		T v0123 = v01 * (1.0f - yf) + v23 * yf;
		T v4567 = v45 * (1.0f - yf) + v67 * yf;
		return v0123 * (1.0f - zf) + v4567 * zf;
	}

	inline T sample(Vector3 v)
	{
		return sample(v.x, v.y, v.z);
	}
};

template<>
void VolumeData<uint8>::updateVolumeTexture()
{
	vtex = Texture::fromMemory("", &data[0], ImageFormat::L8(), dim.x, dim.y, dim.z, ImageFormat::L8(), Texture::DIM_3D_NPOT, Texture::Settings::video());
}

template<>
void VolumeData<uint16>::updateVolumeTexture()
{
	vtex = Texture::fromMemory("", &data[0], ImageFormat::L16(), dim.x, dim.y, dim.z, ImageFormat::L16(), Texture::DIM_3D_NPOT, Texture::Settings::video());
}

template<>
void VolumeData<float>::updateVolumeTexture()
{
	vtex = Texture::fromMemory("", &data[0], ImageFormat::L32F(), dim.x, dim.y, dim.z, ImageFormat::L32F(), Texture::DIM_3D_NPOT, Texture::Settings::video());
}

template<>
void VolumeData<double>::updateVolumeTexture()
{
	char *p = (char *)&data[0];
	for(int i = 0; i < data.size(); ++i)
	{
		int index = i * sizeof(double);
		p[index] = (p[index + 0] ^= p[index + 7]) ^ (p[index + 7] ^= p[index + 0]);
		p[index] = (p[index + 1] ^= p[index + 6]) ^ (p[index + 6] ^= p[index + 1]);
		p[index] = (p[index + 2] ^= p[index + 5]) ^ (p[index + 5] ^= p[index + 2]);
		p[index] = (p[index + 3] ^= p[index + 4]) ^ (p[index + 4] ^= p[index + 3]);
	}

	Array<float> fdata;
	fdata.resize(data.size());
	for(int i = 0; i < data.size(); ++i)
		fdata[i] = float(data[i]);
	vtex = Texture::fromMemory("", &fdata[0], ImageFormat::L32F(), dim.x, dim.y, dim.z, ImageFormat::L32F(), Texture::DIM_3D_NPOT, Texture::Settings::video());
}

#endif
